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ABSTRACT 

A new, much improved model of the Galactic Magnetic Field (GMF) is presented. We use the 
WMAP7 Galactic Synchrotron Emission map and more than forty thousand extragalactic rotation 
measures to constrain the parameters of the GMF model, which is substantially generalized compared 
to earlier work to now include an out-of-plane component (as suggested by observations of external 
galaxies) and striated-random fields (motivated by theoretical considerations). The new model pro- 
vides a greatly improved fit to observations. Consistent with our earlier analyses, the best-fit model 
has a disk field and an extended halo field. Our new analysis reveals the presence of a large, out-of- 
plane component of the GMF; as a result, the polarized synchrotron emission of our Galaxy seen by 
an edge-on observer is predicted to look intriguingly similar to what has been observed in external 
edge-on galaxies. We find evidence that the cosmic ray electron density is significantly larger than 
given by GALPROP, or else that there is a widespread striated component to the GMF. 



1. INTRODUCTION 

Magnetic fields are ubiquitous in the Galaxy. They 
permeate the interstellar medium and extend beyond the 
Galactic disk, and they are present in stars, supernova 
remnants, pulsars and interstellar clouds. The magnetic 
field in the diffuse interstellar medium has a large-scale 
regular component as well as a small-scale random part, 
both having a strength of order micro- Gauss. 

The large-scale Galactic magnetic field (GMF) has re- 
ceived considerable attention yet it remains poorly un- 
derstood. The main difficulty in determining the large- 
scale GMF is the lack of in situ measurements of the 
magnetic field. The best available constraints are Fara- 
day rotation measures (RM) and polarized synchrotron 
radiation (PI), both of which are line-of-sight integrated 
quantities. The RM (PI) depends on the component 
of the field parallel (perpendicular) to the line-of-sight, 
weighted by the total (relativistic) electron density n e 
(n cre ). This complementarity in the sensitivity to or- 
thogonal magnetic field components and different elec- 
tron distributions is a powerful reason for combining the 
two data sets in a joint analysis. 

O ur previous syst e matic effort to combine these data 
sets, ||Jansson et alT] (|2QQ9| hereafter JFWE09), investi- 
gated the validity of the large-scale Galactic magnetic 
field models in the literature at that time, by testing their 
predictions for polarized synchrotron and extragalactic 
rotation measure data. It was found that all extant mod- 
els failed to provide a good fit to the measured RMs and 
PI maps, even when their parameters were re-optimized 
to fit the data: their functional forms were simply not 
general enough to reproduce important features of the 
data. Some simple modifications to existing models were 
investigated in JFWE09 which improved the fit. In par- 
ticular, the magnetic field in the halo was found to have 
a form which is fundamentally different than the field in 
the disk, rather that being a weaker version of the disk 
field. 

In this paper we make use of vastly more RM data 



than previously available, and we update to the latest 
WMAP7 synchrotron emission data. Even more impor- 
tant are the changes we have made to the form of the 
GMF model considered. We allow here for the possibil- 
ity of a large-scale out-of-the-plane component and al- 
low for striated and fully random fields. The need for 
an out-of-the-plane component to the field is suggested 
by ob servations of external galaxies ( |Beck|[2QTT| |Krause| 
2009]). 

The fit to the new, more general field model confirms 
the need for these new components, and the resulting 
GMF gives a dramatic improvement in the quality of the 
fit to the data, even as the quality and quantity of data 
have improved. The results presented here substantially 
revise our understanding of the Milky Way's magnetic 
field. 



Som e notable recent works include Jaffe et al. (2010 



2011 ) studying the Galactic disk field with synchrotron 
data and allowing for "ordered random" magnetic fields; 

in et al-| ( |2008b and JSun fc Reich| ( |2010| ) modeling 
the disk and halo GMF and constraining the model 
with mult i- wavelength synchr otron and rotation measure 
data; |Pshirkov et al. ( |2011| ) comparing models in the 
literature (and proposing two benchmark models) using 
full-sky rotation measure data, some of which is unpub- 
lished. 



2. METHOD 



We use t he numerical Hammurabi code (Waelkens 



et al.||2009[ ) to calculate simulated data sets of Rota- 
tion Measures and the Stokes parameters Q and U, 
from 3D models of n e , n cre and B. As an estimator 
of the quality-of-fit to the parameters of the large-scale 
GMF, we use Xtot = %mXrm + w Qu(Xq + Xu)> wnere 
the coefficient factors iurm,qu are chosen to give equal 
weight to the RM and synchrotron data sets, and, e.g., 
Xq = Ei(Qdata,i - Qmodei,i) 2 /crQ 5 i, where the sum runs 
over the individual pixels. With x?ot a function of GMF 
parameters, we use a Metropolis Markov Chain Monte 
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Carlo (MCMC) algorithm ( [Metropolis et al |1953| ) to find 
best-fit parameters and confidence levels tor the GMF 
model. 



The variances in the observables - Gq • , crfj 



and a\ M i 

- which are needed to evaluate \ 2 are n °t merely the ob- 
servational or experimental errors, but include and are 
dominated by the astrophysical variance caused by tur- 
bulent magnetic fields and inhomogeneities in the inter- 
stellar medium. The estimation of these variances is cen- 
tral to our analysis, and is discussed in the next section. 

3. DATA 
3.1. Faraday rotation measures 
The rotation measure, in units of rad m -2 , is 



RM ~ 



Jo 



n e (l) \ f B\\(l) 



(1) 



where n e is the total density of ionized electrons, which 
is dominated by the thermal electron density. Rotation 
measure is inferred from the relation between the polar- 
ization angle of a source and the wavelength- squared of 
the observation: = Oq + RM A 2 , in the Faraday-thin 
case. The reliability of the estimated RM thus depends 
on the number of and spacing between the wavelengths 
with which the source has been observed. 

The publicly available extragalactic RM data has in- 
creased by more than an order of magnitude since 
JFWE09, t hanks to the r e-an alysis of NVSS polariza- 



tion data by Taylor et al. 



( |2QQ9| ). This data set includes 
37543 RMs that cover the sky north of declination —40°. 
However only two wavelengths were used in the deriva- 
tion of these RMs, so they are the least reliable RMs in 
our sample. Complementing these RMs, we include in 



our analysis 194 recently obtained disk RMs by | Van Eck 



et al. ( 2Q11|); 380 RMs from the Canadian Galactic Plane 



SunFey prown et al.||2Q Q3 |); 148 RMs from the Southern 
Galactic Pla ne Survey (Brown et al. 2007k; 813 high lat- 
itude RMs (|Mao et al. 2010) 60 R Ms near the Small 
Magellanic Cloud (Mao et al.||2008| and 200 RMs near 
the L arge Magellanic Cloud ([Gaensler et al. 2005; [Maol 
|2012| ); 160 RMs near Centaurus A (iFeain et al ||200^ ; 
and 905 RMs from various other observational efforts 



( Simard-Normandin et al. 1981; 


Broten et al.|1988; Clegg 


et al.|iO^||Oren k Wolie|l995; 


Minter & Spangler|1996| 


Gaensler et al.||200ip The total 


number of extragalactic 



To avoid skewing the mean and variance of the rota- 
tion measure for a particular direction, we need to re- 
move data points that are in fact multiple measurements 
of the same s ource. We do this by mapping the RMs to a 
HEALPhf^ Gorski et al.|2Q05| ) pixelation of the sky, with 



> x 10 -4 square-degree pixels (i.e., about 50 million pixels 
for the full sky). Multiple measurements within a single 
pixel are averaged. The various RMs in our combined 
sample have been determined from different numbers of 
wavelength measurements, and can divided into three 



groups of increasing reliability; the |Taylor et al. ([2009 ) 



data is derived using only two wavelengths, Br oten et al.| 
(1988) used a few widely spaced wavelengths, and the 
other RMs in our sample used several closely spaced 

1 http://healpix.jpl.nasa.gov 
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Figure 1. Top: The RM sky, after removing outliers and aver- 
aging to 4 °-bv-4° pixels. Middle: The nearby HI bubble seen in 
RM, from |Wolleben et al.| ( |2010| >. Bottom: The RM sky with the 
nearby RM bubble subtracted. 

wavelengths. When a pixel has multiple RMs from a 
mix of these three groups, only those from the most reli- 
able group are kept, and then averaged. This procedure 
leaves 38627 pixels. 

Plausible outliers - sources with large RM contribu- 
tions likely due to effects other than the Galactic mag- 
netic field (e.g., RM intrinsic to its source) - are removed 
by an iterative scheme: i) For each pixel the mean and 
variance of the RM of neighboring pixels (those within 
2°) are calculated, ii) If the RM of the pixel deviates 
from this mean by more than three standard deviations, 
it is removed. This process is repeated until no RMs are 
marked for removal. In our sample, three such cycles are 
necessary, and results in the removal of a total of 666 
pixels. 

Obtaining an accurate estimate of the astrophysical 
variance due to random magnetic fields and inhomo- 
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Figure 2. Skymaps of observables and fits in Mollweide projection. Galactic longitude / = 0° in the center and increasing to the left. 
Columns, from left: rotation measures (in rad/ra 2 ), Stokes Q and Stokes U (in mK). Rows, from top: data, simulated data from best-fit 
model, cr and the contribution of the pixel to x 2 . White pixels correspond to either missing data (RM), or masked data (PI). The simulated 
RM map also includes predictions for regions without data. 



geneities in the magnetized ISM is crucial to the entire 
analysis. By simulating sky-maps of rotation measures 
using large-scale magnetic field models, such as the one 
used in this paper, we find the rotation measure varies 
only slightly on small angular scales (« few degrees). 
Hence we bin the 37961 pixels to a set of 2670 approxi- 
mately 4°-by-4° pixels. (The full sky has 3072 pixels but 
some portions of the sky have no measured RMs.) The 
sub-pixels contained in each of these larger pixels are 
used to calculate the variance of the rotation measure in 
each large pixel. 

In a few cases, the number of sub-pixels with measured 
RM in a large pixel is less than iV min = 10. In this 
case we successively increase the search radius centered 
on the given pixel up to r = 4°, until iV m i n RMs are 



found. For a small number of pixels, N < iV min even 
when r = 4° , in which cases we de- weight these pixels by 
increasing their estimated variance by a factor N m [ n /N. 
If the variance of any of these pixels is less than the 
average for that meridian, it is replaced by the average. 
This is required for only 12 pixels. If no sub-pixels are 
found within r degrees, that pixel is excluded. We are 
left with 2637 RM pixels in the end. We note that the 
observational uncertainty in the RMs is not explicitly 
included in the calculation of the total variance, since it 
enters implicitly in the variation in sub-pixel RMs and 
moreover the measurement error is small compared to 
the natural variance. 

3.1.1. Foreground subtraction 
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In Wolleben et al. ( |2Q1Q ) the authors do a rotation 
measure synthesis analysis, using the first results from 
the Global Magneto-Ionic Medium Survey (GMIMS), 
and find that the RM of a significant portion (about 
1/20) of the sky is dominated by a local HI bubble. The 
nearby RM contribution fro m the bubble can be se en in 
figure |T] which is taken from Wolleben et al.| ( |2Q1Q[ ) and 
smoothed to 4°-by-4° pixels. The shown region is a circle 
of radius 30° centered on (I, b) = (40°, 30°). 

In our main analysis we perform the model optimiza- 
tion after subtracting the nearby H I bubble from the RM 
data set. However, we also do the model fitting a sec- 
ond time, using the unsubtracted RM data. Performing 
the fit twice provides an important sanity check about 
the Galactic magnetic field model: if the optimized x 2 
is worse when the nearby feature has been subtracted, 
the field model would probably be a poor appro ximation 
to the large-scale GMF. As reported in §7.1.1} the opti- 
mized x 2 is gratifyingly lower when the local H I bubble 
is removed. 

Several other nearby structures exist, and when RM 
synthesis is available for them, their contributions can 
also be subtracted before fitting the global GMF model 
to RM data. We expect GMIMS and surveys similar 
to it, to yield an increasingly accurate map of RM fore- 
grounds. 

3.1.2. Pulsars 

Pulsar rotation measure data should in principle pro- 
vide significant additional constraints on the Galactic 
magnetic field. However, the majority of pulsars have 
poorly estimated distances, so the predicted RMs are 
correspondingly very uncertain, particularly for lines-of- 
sight where the magnetic field has reversals. Properly 
estimating a for pulsar RMs is also less straightforward 
than for extragalactic RMs. Thus pulsars are not used 
in the present analysis but will be included at a future 
stage. 

3.2. Polarized synchrotron emission 

The polarized radiation at 22 GHz is dominated by 
Galactic synchrotron emission. For a power-law distri- 
bution of relativistic electrons (n cre ) with spectral index 
s, the synchrotron emissivity is 



1 + S j 



(2) 

For a regular magnetic field and a power-law distribution 
of electrons with spectral index 8 = 3, the emitted syn- 
chrotron radiation is linearly polarized to around 75%. 
Observationally, the polarization fraction is much lower 
due to depolarizing effects, such as the presence of tur- 
bulent or otherwise irregular magnetic fields which de- 
polarize the radiation through line-of-sight averaging. In 
this paper we will use the polarized components of the 
synchrotron data (the Stokes Q and U parameters) to 
constrain the large-scale magnetic field model. 

In the WM AP seven-year r elease of their K-band (22 
GHz) data QGold et aL 2011) the data is separated into 
foreground components, including synchrotron emission. 
At this frequency we can assume that the Faraday ro- 
tation of the synchrotron data is negligible; thus Q and 
U are independent of RM. We take the WMAP7 syn- 
chrotron data set and average the Stokes Q and U pa- 




Figure 3. T he polarized sync hrotron masks. Top: The black re- 
gion show the |Gold et al.||2011| ) mask (covering 27% of the sky); 
the gray region shows the expanded mask used in the main anal- 
ysis of this paper (covering 35% of the sky). Bottom: Two very 
differ ent m asks derived from the "pull" of the polarized intensity 
(see 33.2.1) ; the black region shows the mask for pull > 3 (cover- 
ing l47o ot the sky); the gray region shows the mask for pull > 2 
(covering 29% of the sky). 

rameters to form HEALPix maps with 4°-by-4° pixels, 
as done for the RM data. The variances of these individ- 
ual pixels are calculated from the original l°-by-l° pixels 
(the resolution of WMAP at 22 GHz). Figure [2] shows 
the processed Stokes parameters. 

3.2.1. Polarization mask 

Synchrotron flux from nearby structures such as super- 
nova remnants pollute the emission caused by the large- 
scale GMF in the diffuse interstellar medium. These 
structures are prevalent in the disk and are best masked 
out in an analysis of the large-scale magnetic field. We 



use the WMAP polarization mask discussed in Gold et al 
(2011 ), covering 27% of the sky, but expand the mask by 
hand to include some partially masked structures and 
remove some distinct high-PI regions that appear to cor- 
respond to localized structures. The final mask covers 
35% of the sky and is shown in Figure [3| We take the 
expanded mask as our primary mask. However, to check 
the sensitivity of our best-fit parameters to the choice of 
polarization mask we also consider the WMAP mask and 
two drastically different masks derived from the "pull" of 
the polarized synchrotron data. For each pixel, we define 

the pull to be p = (Q 2 + /7 2 )/(<7q + <r^), and create 

two masks, for p > 2 and p > 3, respectively. Mask- 
ing out regions with a large p should remove the most 
prominent local structures, such as the Northern Spur. 
However, it may also remove important regions of signif- 
icant PI caused by the large-scale GMF. As can be seen 
in Figure [3] the masks are very dissimilar to our primary 
mask, and this is the main reason we include them. In 
section 7.1 we will see that the conclusions of our analysis 
are not sensitive to the choice of synchrotron mask. 
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4. ELECTRON DENSITIES 

The rotation measures and synchrotron emission are 
line-of-sight integrals of the magnetic field but weighted 
by the thermal and relativist ic (also known as cosmic 
ray) electron densities, n e and n cre , respectively. In this 
paper we adopt t he standard NE2QQ1 the rmal electron 
density model by Cordes & Lazio (2002) for n e , with 
the mid-plan e density and vertical scale-height modified 
according to |Gaensler et al.| ( |2008| ) . 

4.1. Relativistic electron density 

We consider two distinct models for the spatial distri- 
bution of re lativistic elec trons: the one obtained from 
GALPR OP (|Strong et al |2009|) , and the one adop ted by 
WM AP ([Page et al.| (I20U71T" who were following Drim- 
mel fc Sp ergel ( 2001 ) J~ The models are fundamentally 
different in that the WMAP model is just a simple phe- 
nomenological parameterization while the GALPROP 
distribution is based on the distribution of supernovae 
remnants in the Galaxy and numerical simulation. The 
GALPROP model is not peaked at the Galactic center 
and is not described by a simple function; we thank A. 
Strong for providing it to us as a FITS file. Both models 
are shown in Figure [4j 

The WMAP modelis 

C cre (r,z) = C cre , e- r//v sech 2 (z/h z ). (3) 
The quantity C cre (r, z) is defined by 



iV(7, r, z)dj = C cre (r, z)j p dj, 



(4) 



where N is the number density. The normalization factor 
C cre ,o is such that for 10 GeV electrons, C cre (Earth) = 



trons at Earth (Strong et al. 2007). We consider two 
variants on the WMAP model: Erst, using the original 
WMAP parameter values, h r = 5 kpc and h z = 1 kpc, 
and second, allowing h r and h z to be free parameters to 
be varied along with the parameters of the GMF in the 
parameter optimization. 

For all models, the number density for other energies is 
calculated assuming a power law distribut ion with spec- 
tral index p — —3 (Bennett et al. 2003). The spatial 
distributions of these models are shown m Figure [4j 

5. GALACTIC MAGNETIC FIELD MODEL 

The most familiar components of the Galactic mag- 
netic field are the large-scale regular fields and the small- 
scale random fields. The latter are due to a vari- 
ety of phenomena including supernovae and other out- 
flows, possibly compounded by hydrodynamic turbu- 
lence, which are expected to result in randomly-oriented 
fields with a coherence leng t h A of order 100 pc or less 



Haverkorn et al 



2008). In 



( jGaensler" fc Johnston||1995| 
addition to these, we include m our model "striated ran- 
dom fields - fields whose orientation is aligned along some 
particular axis over a larger scale, but whose strength 
and sign varies on a small scale. Such striated fields can 
be produced by the levitation of bubbles of hot plasma 
carrying trapped randomly oriented fields away from the 
disk, or by differential rotation of small scale random 
fields, or both. The predominant orientation of stri- 
ated fields produced by differential rotation is plausibly 
aligned with the local coherent field. Striated fields are a 



o 

CL 




Figure 4. Top: The spatial distribution of relativistic electrons 
used in our main analysis (from GALPROP | Strong et al.| ( |2004| 
|2010| and private communication)). The contour units are cm"" 3 
tor 10 GeV electrons. Middle: GALPROP n nC re increased by a 
factor of 2.9, which optimize our % 2 under the a ssum ption that 
no striated fields are present in the Galaxy (see j |6.2| ). Bottom: 
Original WMAP n nC re, with radial and vertical scale height, 5 kpc 
and 1 kpc, respectively. 

special case of the more gen eric possibility of an isotropic 
random fields introduced in Sokoloff et al. (1998), which 
can be considered a superposition ol multiple striated 
and purely random fields. 

These three distinct types of magnetic structures - 
large-scale regular fields, striated fields, and small-scale 
random fields - can be disentangled because they con- 
tribute differently to different observables. The large- 
scale regular field contributes to all the observables - I, 
PI and RM - while the small-scale random field only 
contributes to the total synchrotron emission, I. In the 
present work, we restrict our analysis to striated and reg- 
ular fields and therefore do not fit I or include the small 
scale random fields in our model. 

The striated field contributes to I and PI, but in leading 
order it does not contribute to rotation measures due to 
its changing sign, except possibly for a very small num- 
ber of pixels for which the line-of-sight is precisely aligned 
with the direction of the striated field, since we smooth 
over pixels whose size is large compared to the coher 
ence le ngth characterizing the field reversals. 
(2010) use the term "ordered random fields" 



( Jaffe et al. 
or what is 

probably phenomenologically equivalent to our striated 
fields - they define it as a field component contributing to 
I and PI but not RM - although their cartoon indicates 
the coherence length for reversals is similar in all direc- 
tions whereas we envisage an origin which would natu- 
rally lead to asymmetric coherence lengths. Since fields 
can be random in some respects and ordered in other re- 
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spects, in a variety of ways, e.g., coherence length could 
depend on direction but field orientations be random, 
we prefer the more vivid and specific term "striated" to 
the term "ordered random", for the type of field being 
described here.) 

5.1. Large-scale regular field 

The necessity of separate disk and halo fields was 
show n in JFWE Q 9, and obser vations of external galax- 
ies (Beck (2009), Krause ( 2009[ )) prompt the inclusion 
of an out-ol-plane held component. Thus we model the 
large-scale regular GMF with three separate components. 
Furthermore, we restrict ourselves to functional forms 
such that each component of the field is separately di- 
vergenceless so their parameters can be specified inde- 
pendently. Imposing flux conservation has not been uni- 
versally adopted in past modeling, because the constraint 
is so restrictive: it can be difficult to find phenomenologi- 
cally appropriate forms which can be explicitly expressed 
in closed form and which are manifestly divergenceless. 
However flux conservation is an extremely important and 
constraining theoretical condition, so we demand that it 
be enforced. 

We use right-handed Cartesian (x, z) and cylindri- 
cal (r, </>, z) coordinate systems throughout the following 
discussion, where the Galactic center is at the origin, 
Galactic north is in the positive z-direction, and the Sun 
is located at x = —8.5 kpc. The field is set to zero for 
r > 20 kpc and in a 1 kpc radius sphere centered on the 
Galactic center. 

5.1.1. Disk component 



For the d isk, we use a generalized form of the Brown 



et al. (2007) model, which is partially based on the struc- 
ture of the NE2001 thermal electron density model. The 
main focus of the present work is on the halo field, so 
we satisfy ourselves with adopting this pre-existing form, 
but adjusting the field strength parameters and depen- 
dence on z to enforce flux conservation and improve the 
fit. 

The disk field is constrained to the a>?/-plane, and de- 
fined for Galactic radii r between 3 kpc and 20 kpc. 
In the 'molecular ring', between 3 kpc and 5 kpc, the 
field is purely azimuthal with a field strength of b r [ ng . 
Between radii 5 kpc and 20 kpc there are eight log- 
arithmic spiral regions with opening angle i = 11.5°. 
The dividing lines between these spiral regions follow 
the equation r = r- x exp(0tan(9O° — z)), where r- x = 
5.1, 6.3, 7.1, 8.3, 9.8, 11.4, 12.7, 15.5 kpc are the radii 
where the spirals cross the negative x-axis. The mag- 
netic field direction in the spiral regions is given by 
b = sin(z)f + cos(i)0. The field strength, ^, in mag- 
netic spiral i is defined at r = 5 kpc, and falls off as 
r -1 . To conserve magnetic flux the field strengths of 
seven of the spirals are free parameters in the model, 
with the field strength of the last spiral set by the 

constraint b% = — ^I=i/A//8> where fa is the rela- 
tive cross-sectional areas of the spirals (for a fixed ra- 
dius). From r_ x above, we can derive the corresponding 
fi = 0.130, 0.165, 0.094, 0.122, 0.13, 0.118, 0.084, 0.156. 

The extent of the disk field is symmetrical with respect 
to the mid-plane, and set by the height parameter ftdisk? 
where the disk field transitions to the toroidal halo field. 



The transition is given by the logistic function, 

L(z,h,w) = (l + e- 2 ^l-^/™)~\ (5) 

where the free parameter Wdisk sets the width of the tran- 
sition region; for small w, L becomes a step function. The 
disk component is multiplied by (1 — L(z, ftdisk? ^disk)) 
and the halo field is multiplied by L(z, ftdisk? ^disk)- 

5.1.2. Toroidal halo component 

The halo field has a purely toroidal, i.e. azimuthal, 
component defined as 

Bf r (r, z) = e-^ z °L(z, /i disk , w disk ) (6) 

'B n (l - L(r, r n , iu h )), if z > 
K B a (l-L(r,r a ,w h )), if z < 0. 

This halo field has an exponential scale height, and sep- 
arate field amplitudes in the north and south, B n and 
£? s , respectively. The northern (southern) radial extent 
of the halo field is set by r n (r s ). The parameter con- 
trols the width of the region where the halo field is cut 
off. 

We considered several forms for the halo field, includ- 
ing axisymmetric and bisymmetric spirals, and settled on 
the purely toroidal model when it was clear that it led to 
a superior fit to data. Some alternative halo co mpo nents 
that we tested, and rejected, are discussed in § A. 1 



5.1.3. Out- of -plane component 

The halo field is generalized compared to earlier work, 
by including an out-of-plane component. We refer below 
to the out-of-plane halo component as the "X-field" com- 
ponent, since it is partially motivated by the X-shaped 
field structures s een in radio observation s of external, 
edge-on galaxies ( |Krause(2009| |Beck||2QQ9| ) . 

We choose the out-of-plane component to be axisym- 
metric and poloidal, i.e., lacking any azimuthal compo- 
nent (which is incorporated via the toroidal halo com- 
ponent). To find a reasonable functional form for such 
a field, that is also divergenceless, is not simple. We 
developed the parametrization below; a visualization is 
provided in Figure [6] for the parameters of the best-fit 
GMF (see Table [lj. The field at any position (r,z) is 
specified, as discussed below, in terms of r p , the radius 
at which the field line passing through (r, z) crosses the 
mid-plane (z = 0). 

We take the field outside galactocentric radius r x to 
have a constant elevation angle, © x , with respect to the 
mid-plane. Within this radius, the elevation angle Ox 
is linear in the radius, becoming vertical, Ox = 90° , at 
r = 0. We define the field strength in the mid-plane by 
the function 

&x(r P ) = Bxe-^ /rx , (7) 

where Bx is the overall amplitude of the X-field and r p is 
the mid-plane radius of the field line that passes through 
(r,z). 

With this general geometry, the requirement V • B = 
is sufficient to fully characterize the field. The field 
line with r p = r x marks the border between the region 
with constant elevation angle and the interior region with 



varying elevation. In the constant elevation region, the 
field strength is &x(fp) ^pA*> where 



M/tan(e°). 



(8) 



In the region with varying elevation angle the field 
strength is instead bx(r p )(r p /r) 2 , and the elevation angle 
and r p are given by 



r = 

p r£ + |z|/tan(e°)' 

W 



Qx(r, z) =tan 1 ( - 
V' 



(9) 
(10) 



Altogether, the out-of-plane component has 4 free pa- 
rameters: £?x, ©x> r x anc ^ r x- 

5.2. Striated random fields 

We include the possibility of striated magnetic fields 
by adding a multiplicative factor to the calculation of 
PI, such that when this factor is equal to unity the model 
describes a purely regular field. We parametrize striated 
and purely random fields as B^ tri = j3B^ eg . We let the 
factor be a free parameter in the large-scale GMF model. 
We originally performed the analysis allowing the disk, 
toroidal halo, and X-field each to have a separate amount 
of striation (see appendix |A| . We did not find a signifi- 
cant improvement in x 2 using this added freedom, so for 
the final parameter optimization used a single /3 value 
for all components. This means the striated field is ev- 
erywhere aligned with the local large-scale field and has 
the same relative magnitude everywhere in the Galaxy. 

When the striated field is aligned with the regular field, 
there is an obvious degeneracy between the strength of 
the striated magnetic field component and the relativis- 
ts electron density: if we write the multiplicative fac- 
tor as 7 = a(l + /?), we can interpret a as being a 
rescaling factor for the relativistic electron density, with 
£>stri = f3B^ eg . The distribution of relativistic electrons 
in the Galaxy is not well enough known to permit this de- 
generacy to be disentangled at present. Of course, since 
f3 > it follows if 7 is found to be less than unity we can 
conclude that a < 1, and that n cre has been underesti- 
mated. 

5.3. Parameter Estimation 

As noted in JFWE09, avoiding false x 2 minima when 
optimizing a model is very difficult, and we have devoted 
considerable effort to exploring the very large parame- 
ter space available for the model outlined in the previ- 
ous section. The model optimizatio n is done using the 
PyMC package by Patil et aT| ( 2Q1Q[) , and uses an adap- 
tive Metropolis MCMC algorithm. To achieve good mix- 
ing and convergence of the Markov chain, we continue 
to sample the parameter space until t he Gelman-Rubin 
convergence and mixing statistic, R ( [Gelman fc Rubin] 

1992 ), satisfies the condition R < 1.03 for all parame- 
ters. 



The final Markov chain has 100k steps, and the 
Monte Carlo standard error for any given optimized pa- 
rameter is at least an order of magnitude less than the 
estimated confidence range of the same parameter. 

6. RESULTS 
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Figure 5. Top view of slices in the x-y-plane of the GMF model. 
Top row, from left, slices at z = 10 pc and z = — 10 pc. Bot- 
tom row, slices at z = 1 kpc and z = — 1 kpc, respectively. The 
color scheme shows the magnitude of the total regular field, with 
negative values if the azimuthal component is oriented clockwise. 
The location of the Sun at x = —8.5 kpc is marked with a circle. 
From the top panels it is clear that the magnetic field just above 
and below the mid-plane are very similar, but not identical, due 
to the superposition of the z-symmetric disk field component with 
the ^-asymmetric toroidal halo component. At \z\ = 1 kpc the field 
is dominated by the halo component, but still exhibits signs of the 
superposition with the X-field, and even the disk field. 




5 kpc 



12 3 

Figure 6. An x — z slice of the galaxy showing only the out-of- 
plane "X" component. The black lines crossing the mid-plane at 
±4.8 kpc traces the boundary between the outer region with con- 
stant elevation angle, and the inner region with varying elevation 
angle. The black arrows show the direction of the field. 

6.1. Optimized large-scale magnetic field model 

The large-scale Galactic magnetic field model has 21 
free parameters. Table[l]lists the best-fit values and 1 — a 
confidence intervals. 

6.1.1. The disk field 

The best-fit field in the disk is shown in the top panel 
of Figure [5] The innermost arrow refers to the molecular 
ring region; consecutive arrows are positioned in spiral 
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Table 1 

Best-fit GMF parameters with 1 



- a intervals. 



Field 


Best fit Parameters 


Description 


Disk 


0\ = U.l ± 1.8 IIKj 


field strengths at r = 5 kpc 




02 = o.U ± U.O flLx 






03 = —0.9 ± 0.8 /iG 






04 = —0.8 ± O.J /iG 






05 = —2.0 ± 0.1 fiKj 






06 = — 4.Z ± 0.5 flij 






67 = 0.0 ± 1.8 /iG 






Og = Z. ( ± 1.8 /it^r 


mierrea trom 01, 07 




°ring — 0.1 ± 0.1 /iG 


ring at 3 kpc < r < 5 kpc 




Aldisk — 0.40 ± O.Oo kpc 


disk/halo transition 


— — — — - — 


^disk — U.Zl ± U.U8 KpC 


transition width 


Toroidal 


-D n — 1.4 zt U.l fJAj 


northern halo 




R ——114-01 //H 


\1 1 f- r 1 OT' T 1 V| 1 

oUUlllclIl IldlL) 




r n = 9.22 ± 0.08 kpc 


transition radius, north 




r s > 16.7 kpc 


transition radius, south 




w h = 0.20 zb 0.12 kpc 


transition width 




zo = 5.3 zb 1.6 kpc 


vertical scale height 


X halo 


B x =4.6zb0.3/iG 


field strength at origin 




e° = 49 zb i° 


elev. angle at z = 0, r > r x 




r x = 4.8 zb 0.2 kpc 


radius where Ox = B x 




r x = 2.9 zb 0.1 kpc 


exponential scale length 


striation 


7 = 2.92 zb 0.14 


striation and/or n cve rescaling 



Note. — For the parameter r s only a lower 68%-bound is given. 
The Markov chain parameter distribution for this parameter, and 
a few others of interest, are shown in Figure [T] 




(kpc) 



w h (kpc) 
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Figure 7. MCMC histograms for a selection of the GMF pa- 
rameters. The counts are in the units of 10 3 . The top left panel 
(frring) shows a Gaussian-like distribution, and is typical for most 
of the parameters in the fit. The cases with significant deviations 
from Gaussianity are shown in panels 2-4. The scale height of the 
toroidal halo component, zo, is close to gaussian, but has positive 
skew. We note that r s is unconstrained for large values. 

arm regions 1 to 8. Because of the superposition of the 
disk with the toroidal halo and X-field, parts of the field 
in the disk become asymmetric in z (e.g, arm region 1 
and the molecular ring). The smooth transition between 
the disk and halo fields is centered around 400 pc, but 
the transition width is large enough that the total field 
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■ regular (x = -8.5 kpc) 

■ regular + max striated (x = -8.5 kpc) 
- regular (x = -10 kpc) 

■ regular + max striated (x = -1 kpc) 
rms B from synchrotron (Cox 2005) 
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Figure 8. The field strength of the optimized GMF model as a 
function of z, at (x,y) = (—8.5,0) kpc (the Solar neighborhood) 
and at (x,y) = ( — 10,0) kpc. The solid lines show the magnitude 
of the regular field and the dashed lines show the magnitude of 
the combined striated and regular field. The dotted line shows 
an es timate of th e total field (including small-scale random fields) 
from |Cox| ( |2005| ). The large difference in predicted field strength 
for x = —W.b kpc and x = — 10 kpc at small \z\ is due to the points 
being located in two different magnetic spiral arms. 

is a mixture of both, even at the mid-pl ane. 

In agreement with Brown et al. (2007) we find a large- 
scale reversal between the Scut urn- Crux spiral arm (re- 
gion 2; counterclockwise field as viewed from the north 
Galactic pole) and the Perseus s piral arm (region 6 ; 
clockwise field). In contrast with Brown et al.| 



we find evidence of another reversal between the Perseus 
and Norma spiral arms (regions 6-8). However, the field 
strength in arm region 8 is less than two standard de- 
viations from zero, hen ce the evidence for t his reversal 
is weak. We note that |Brown et al. ( |2007[ ) only mod- 
eled the GMF in the region 253 u < / <3H8 U , and would 
thus not be very sensitive to data constraining region 7 
and 8. We also note that |Brown et al.| ( |2007[ ) reported 
a counterclockwise field in the molecular ring, while we 
find a very weak field that is mostly present in t he model 
via the supe rimposed halo a nd X-field. In the |Van Eck 
et al. ( 2011[ ) extension of the |Brown et al. ( 2007| ) model, 
the authors split the model molecular ring into two half- 
rings, and find a preference for their magnetic fields to go 
in opposite directions. Since this configuration violates 
the divergenceless condition we did not consider such a 
feature in our model. The van Eck fit could be a hint 
that the magnetic field in the molecular ring is not as 
simple as a purely azimuthal field, and explain why in 
our optimized model the field is essentially nonexistent. 

In the past, there has been much discussion of the num- 
ber of field reversals in the disk. In this new model, due 
to parameters having error assignments and there be- 
ing multiple components contributing to the field at any 
given point, the question must be made more precise. 
One could for instance identify loci at which the sign of 
the field differs between adjacent regions in which the 
fields are non-zero, by at least 3cr; simply counting the 
times the arrows in Figure [5] change direction is not suf- 
ficient. 

It is important to stress that the particular functional 
form of the logarithmic spiral has highly non-local impli- 
cations. In reality the observables mostly constrain the 
magnetic field within several kpc of our position, so that 
one should not take too seriously the predictions for the 
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magnetic field strength and orientation on the other side 
of the Galaxy. The disk field deserves further elabora- 
tion and exploration in future work, for instance allowing 
the field to vary smoothly across the arms, allowing the 
arms' locations and widths to be free parameters, and, 
most importantly, exploring alternat ives to the logarith- 
mic spiral structure. For instance, Moss et al. (2012) 
recently showed how dynamo action can lead to the disk 
field on one side of the Galaxy being drastically differ- 
ent to the field on the other side, not having a particu- 
larly orderly form, etc. Devising a way to characterize 
such possibilities, in a way that allows the field struc- 
ture to be constrained, will be a central goal for the next 
phase of this modeling program. Using RMs for pulsars 
with accurately known distances, simultaneously fitting 
the electron densities, and considering PI at different fre- 
quencies, will all have an important role to play in future 
improvements to the disk modeling. 

6.1.2. The toroidal halo field 

A slice through z = ±1 kpc (Figure [5| shows mainly 
the toroidal halo field. In the northern halo, the field 
extends to r ~ 9.2 kpc, while the southern component 
stretches farther, to r > 16 kpc (the actual value is un- 
constrained for large radii). The halo field has a small 
transition width, w 0.2 kpc. Figure [5] also shows the 
significant impact the X-field has on the magnetic field 
in the x-?/-plane by causing the effective pitch angle to 
vary with radius. The toroidal halo field itself has zero 
pitch. 

The most difficult quantity to constrain in the GMF 
model is its vertical extent. In our specific model this is 
mainly set by the scale height of the toroidal halo field. 
This parameter is very sensitive to the chosen electron 
distributions, and indeed for the original WMAP n cre 
only a lower bound on the toroidal scale height is found. 
With n cre from GALPROP the scale height of the field 
is constrained, and is found to be zq = 5.3 ±1.6 kpc. 

Note that because our field model has several compo- 
nents, its vertical profile cannot simply be characterized 
by a single vertical scale-height. This is illustrated in 
Figure |8j which shows the magnitude of the field as a 
function of z for our projected planar position (x = —8.5 
kpc, y = kpc) and a second point 1.5 kpc farther out 
on the x-axis. For comparison, the plot also shows an es- 
timate of the total field magnitude including the random 
com ponent (infe rred from synchrotron emissivity, taken 
from [Cox] ([20051)). 



6.1.3. The out- of -plane field 

The optimized out-of-plane component is significant in 
both strength and extent, and does in fact exhibit an 
"X"-like geometry. The field orientation and strength is 
shown in Figure [6] The field transitions from a constant 
angle to a linearly increasing angle at around 5 kpc. In 
the outer region the elevation angle is approximately 50 
degrees, and the elevation increases at smaller radii until 
the field is completely vertical at r = 0. 

6.2. Striated fields and relativistic electrons 
We optimized t he G MF model with the two n cre mod- 

cre g^eS 



WMAP distribution improves the fit to x 2 /dof = 1.101. 
The WMAP n cre is plotted in the lower panel of Figure 
Hj The GALPROP distribution gives a slightly better 
fit with a reduced x 2 of 1.096. Since the WMAP dis- 
tribution has two free parameters and the GALPROP 
distribution has none, and the GALPROP distribution 
is constrained by a variety of other data, we adopt the 
GALPROP n cre model. 

The best fit value of the product of the striation con- 
tribution and relativistic electron de nsity rescaling is 
7 = a (1 ± /?) ~ 2.9. As noted in £5.2, the present 
analysis does not allow us to discriminate between these 
two sources of increased polarized synchrotron emission. 
A third possibility is that the thermal electron density, 
n e , has been overestimated. In this case, using the cor- 
rect n e would require a stronger GMF to account for 
the observed rotation measures, which in turn would de- 
crease the need for striated fields (or increased n cre ) to ac- 
count for the polarized synchrotron intensity. Of course, 
a combination of all three effects may at work. However, 
since the thermal electron density is a more carefully con- 
strained quantity than the relativistic one, we consider 
it more likely that the large 7 should be interpreted as 
an indication of striated fields in the Galaxy and/or that 
n cre is underestimated. 

Figure [8] shows the contribution of striated fields, if 
the GALPROP and NE2001 models of the electron den- 
sities are correct. The middle panel in Figure [4] shows 
the rescaled GALPROP n cre , under the assumption that 
there are no striated fields in the Galaxy, and the large 7 
is instead due to an underestimated relativistic electron 



els described in §4.1| The original WMAP n t 
a poor fit but optimizing h r and h z appearing in the 



density. Using the parametrization defined in [5.2, the 
n cre is in this case underestimated by a factor a = 7 = 
2.92 ±0.14. 

To investigate further the degeneracy between striated 
fields and the relativistic electron density (and the sen- 
sitivity of our best-fit magnetic field parameters to the 
uncertainty in n cre ) we made the following test: we re- 
optimize the field parameters after multiplying the GAL- 
PROP n cre by a factor exp(|z|/z cre ), with z cre = 10 
kpc. This multiplicative factor increases the effective 
scale height of the relativistic electrons (the number of 
electrons increase approximately by, e.g., 10% at \z\ = 1 
kpc, and by 20% at \z\ =2 kpc). The best-fit parameters 
change on average by 0.4 standard deviations, with most 
of the change predictably being in a, which decrease to 
2.65. The best-fit parameters for the disk field, and geo- 
metric quantities such as r s , r n , ?ib anc ^ ®x are all 
essentially unchanged. The best-fit model is thus robust 
under this degree of uncertainty in n cre . 

This rescaling (with z cre = 10 kpc) also slightly im- 
proves the fit of our model, and could be a sign that 
the GALPROP n cre underestimates the scale-height of 
relativistic electrons. We note that the significant ver- 
tical fields present in our model would tend to increase 
the diffusion of relativistic electrons to larger \z\. GAL- 
PROP currently does not include anisotropic diffusion 
in the calculation of its electron density model, which 
would be necessary to take this effect into account, how- 
ever. Finally, we note that on physical grounds, having 
a striated field which is equally important everywhere 
in the Galaxy (as implied by the best-fit j3 values be- 
ing the same for all three field components) seems some- 
what implausible, favoring the interpretation of a need 
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Figure 9. The Milky Way as seen (in polarization) by an ex- 
tragalactic observer, face-on (above) and edge-on (below). Plot- 
ted "bars" (sometimes referred to as "vectors") are the would- 
be-observed polarization angles, rotated 90° to line up with the 
magnetic field orientation. Lengths of bars are proportional to 
polarization intensity. Faraday depolarization and beam depolar- 
ization are neglected. The face-on plot is overlaid on the NE2001 
thermal electron distribution. 

for rescaling the relativistic electron density rather than 
a large striated field. 

In future work we will incorporate the electron densi- 
ties self-consistently in the overall GMF modeling, and 
generalize GALPROP to include anisotropic and spa- 
tially varying diffusion when calculating n cre . 

6.3. The Milky Way to an external observer 

A hypothetical view of the Milky Way in polarized 
radio emission as seen by an extragalactic observer is 
shown in Figure [9j These maps c an be com pared to ex- 
ternal galaxies presented in, e.g., Beck et al. (2002, see 
also http:/ /www. mpifr-bonn.mpg.de/staff/wsherwood/ 
mag-fields.html for an atlas of magnetic fields in nearby 
galaxies, compiled by R. Beck and W.A. Sherwood). The 
polarization bars (rotated 90° to be aligned with the 
magnetic field direction) are overlaid on the NE2001 elec- 
tron density model. The face-on view shows a tightly 
wound spiral pattern, mostly aligned with the matter 
spiral arms. This outcome was not a foregone conclu- 
sion, since the superposition of the three large-scale field 
components could in theory yield radically different con- 
figurations. 

The edge-on view shows a strong resemblance to the 
polarization patterns seen in some external galaxies, such 
as NGC 891 and NGC 5775, shown in Figure fTOl Mag- 




Figure 10. To y: The magnet ic structure of Milky Way analogue 
NGC 891, from iKrausel 12009), with permission. Contours show 
the total radio intensity, the bars show the magnetic field orien- 
tation (copyright: MPIfR Bonn). The radio map is overlaid on 
an optical image from Canada- France-Hawaii Telescope/ (c) 1999 
CFHT/Coelum. Bottom: The spiral galaxy NGC 5775. Contours 
and bars are again total radio in tensity and magnetic field orien- 
tation. From [Soida et aT] ( |201l] ), with permission. The physical 
width of the field ot view is approximately 30 kpc for both galaxies. 

7.1. Quality- of -fit 

The reduced x 2 of the best fit for the global GMF 
model described above is 1.096, with 6605 data points 
and 21 free parameters. This is a substantial improve- 
ment in fit over previous models, which have reduced 
X 2 > 1-3. 

We note that x 2 serves as a figure-of-merit to compare 
the quality- of- fit for parameter estimation. We have not 
taken steps to assure that the absolute value of x 2 as de- 
fined has the meaning attached in a x 2 -distribution. In 
particular, the low signal-to-noise in parts of the polar- 
ized synchrotron data leads to slightly inflated ctq and 
gjj which we have not corrected because it does not im- 
pact parameter estimation and the ability to compare 
different models' relative fit to the data. 

7.1.1. Sensitivity to foregrounds and choice of synchrotron 

mask 

Performing the parameter optimization withou t first 



netic fields similar to the out-of-plane component de- 
scribed in this paper could thus be present in galaxies 
such as NGC 891 and NGC 5775. 

7. DISCUSSION 



subtracting the nearby H I bubble as discussed in §3.1.11 
leads to a worse x 2 P er degree of freedom: 1.110 instead 
of 1.096 (the best-fit parameters are changed, on aver- 
age, by 1.1 standard deviations). Because only a small 
fraction of the data points are affected, this change in 
the total x 2 is quite significant. Since we should expect 
a correct global GMF model to give a better fit to the 
data if a foreground contaminant is removed, this adds 
credence to our model being correct. 

Performing the optimiza tion with the less conservative 
WMAP polarization mask ( |Gold et al.|2011[ ), the reduced 
X 2 is notably higher, at 1.243. The best-fit parameters 
are quite robust, however; they change on average by 
2 standard deviations. The largest impact occurs for 
the parameter £>x, which changes from 4.6 ± 0.3 /iG to 
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6.9 ±0.4 /iG, which is a significant change in terms of the 
number of standard deviations, but does not substan- 
tially alter the field (e.g., the geometry and extent of the 
X-field does not change by much). 

Optimizing the model with the two masks derived from 
the pull of the polarized intensity (see Figure [3| yields 
best-fit parameters that are on average less than two 
standard deviations from our quoted values. We con- 
clude from this, that while our parameter optimization 
is indeed sensitive to the choice of synchrotron mask, our 
general results for the new model are robust. 

7.1.2. Model comparison 

For model comparison, we opti mize the 1 1 free parame- 
ters of the S un et al. GMF model (ISun et al.|2008||Sun & 



Reich|[2QT0l ), to give the best fit to our observables. The 



optimized Sun et al. model has X? e duced = 1-325, com- 
pared to 1.096 for our model. We note that due to the 
large number of degrees of freedom (6605) the difference 
in X?educed between the two models is truly substantial. 

A figure-of-merit that penalizes for the number of 
model parameters is the Bayesian Information Criterion, 
defined as BIC = Xtot + &l°g(^0> where k is the number 
of free parameters and N is the number of data points. 
For our best-fit model, BIC = 7401, and for the opti- 
mized Sun et al. model BIC = 8832. The major differ- 
ences between our model and the Sun et al. model are 
our inclusion of an out-of-the-plane component and our 
inclusion of striated random fields (or a rescaled n cre ), 
both of which significantly improve the fit. 

As a final compariso n, we also opti mize the bisymmet- 
ric spiral (BSS) field flStanev|[l997l ). The BSS field is 
still often used in the literature, e.g., to predict ultra- 
high energy cosmic r ay deflections ( Takami fc Sato|2010[ 
Vorobiov et al.|2009 among others) , despite having been 



shown to be a poor fit to data ( |Sun et aT||2008| |Jans-| 
|son et aL||2009[ ), and we find the optimized BSS model 
does significantly worse than our model or the Sun et 
al. model, with X? e duced = 1-777 (7 free parameters) and 
BIC = 11790. 

7.2. UHECR deflections 

In our model, ultrahigh energy cosmic rays (UHECRs) 
are deflected predominantly by the toroidal halo field and 
the X-field component, apart from UHECRs observed 
in a direction close to the Galactic plane. Due to the 
asymmetric nature of this field, the average UHE proton 
deflection in the southern part of the Galaxy is approxi- 
mately 60% larger than in the north. For a 60 EeV pro- 
ton, the average deflection (across the sky) is 5.2°, with a 
quarter of the sky having less than 2.2° deflections. The 
magnitude of the deflection is highly non-uniform across 
the sky. 

The predicted deflection for a 60 EeV proton is shown 
in Figure [TTJ Three things of note are apparent from 
the figure: i) the predicted deflection - proportional to 
the integrated transverse magnetic field along the tr ajec- 
tory - differ s greatly bet ween o ur model and that of |Sun| 



et al | fl2008|) and |Stanev d!997[ ). ii) Our predicted denecT 
tion is highly asymmetric across the sky. iii) We predict 
generally larger deflections. The X-field, in particular, 
has a significant impact on the predicted deflections in 
directions towards the inner part of the Galaxy. 





Figure 11. From the top, predi cted deflection an gles for 60 EeV 
protons for our best fit model, the Sun et al. ( 2008 ) model, and the 
bisymmetric spiral model of Stanev ( 199y|. The plots are in the 
Mollweide projection, with Galactic longitude increasing to the left. 
UHECR deflection is proportional to the strength of the magnetic 
field transverse to the UHECR propagation direction, and thus 
provides a useful metric by which to compare different magnetic 
field models. 

Finally, we note that the deflections predicted by the 
best-fit parameters obtain ed u sing the rescaled GAL- 
PROP electron density in §6.2| only differ on average by 
0.3° from the the above case. 

8. SUMMARY AND CONCLUSIONS 

We have developed an improved model for the Galac- 
tic Magnetic Field, whose parameters are determined by 
fitting a large number of Faraday Rotation Measures 
and Polarized Synchrotron emission data. We use the 
WMAP7 maps of synchrotron emission, and rotation 
measures of 40403 extragalactic sources, smoothed on 
4° x 4° pixels, to arrive at 6605 independent observables. 
A key element of our procedure is to determine empiri- 
cally the value of a for each observable, from the variance 
in the observations within each 4° x 4° pixel, giving the 
proper relative significance for each data point in the fit. 

The new 21-parameter GMF model is fundamentally 
different from and more general than any GMF model 
considered previously in the literature. Flux conserva- 
tion is enforced separately for each component and pro- 
vides a powerful implicit constraint, in addition to the ex- 
plicit constraints from fitting the observables. The GMF 
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obtained here gives a far better fit to the observables than 
previous models. Our GMF has a % 2 /dof = 1.1 com- 
pared to y 2 /d of = 1.3 for the best previous form, t he Sun 
et al. model ( [Sun et al.||2QQ8| |Sun fc Reichl|2010| ), after 
optimizing that models' parameters to give trie best fit to 
the same set of observables. The dramatic improvement 
is due to two factors. First, we developed and included a 
closed-form expression for a divergence-free out-of-plane 
field with a sufficiently general, phenomenologically ap- 
propriate geometry. Second, we allowed for the presence 
of large-scale striated fields, or a rescaling of the assumed 
relativistic electron density. 

We adopted for the disk field the general log-spiral form 
that has been used by others, modifying the parameters 
to enforce flux conservati on and optimize the fit. We 
confirm our earlier result ( |Jansson et al.| |2009) that the 
toroidal component of the halo held has its own features 
and cannot be described as a simple scaling of the disk 
field; among other differences, there is an asymmetry be- 
tween the properties of the toroidal component in the 
northern and southern hemispheres, 

We have explored the sensitivity of our results to var- 
ious assumptions, and find that the inferred model pa- 
rameters are generally quite robust. Different choices of 
masks, based on different criteria, do not change the re- 
sultant best-fit GMF models very much. In the present 
work, we adopted the standard Cordes-Lazio NE2001 
thermal el ectron density, with s cale height revised ac- 
cording to Gaensler et al. ( |2QQ8[ ). We considered both 
the GALPROP relativistic electron distribution and also 
the WMAP double-exponential form with two free pa- 
rameters which we fit; the GALPROP distribution is 
physically motivated, has no free parameters and gives a 
better fit, so we adopted it. 

Besides the greatly improved fit, two additional pieces 
of evidence give confidence in the main features of this 
new model and vindicate our methodology: 1) The syn- 
chrotron emission of the Milky Way seen by extragalac- 
tic observers, predicted using the new GMF, resembles 
rather closely observations that have been made of exter- 
nal galaxies from both face-on and edge-on perspectives. 
2) The fit to data improves when the RM of a nearby H I 
bubble is removed, while the best-fit large scale field does 
not change significantly. The fit also prefers the GAL- 
PROP n cre over less physical alternatives considered. 

In future work we plan to do a simultaneous fit, con- 



straining self-consistently parameters of the thermal and 
relativistic electron densities along with the parameters 
of the GMF. Another future direction is to model the 
most important local structures; this should reduce still 
further the small current sensitivity to masks, and pro- 
vide valuable detail about the local ISM. Better knowl- 
edge of the local environment will also benefit the deter- 
mination of global properties from line-of-sight measure- 
ments: due to the fact that all lines-of-sight penetrate 
the local medium, if the local value of n e , n cre or mag- 
netic field is substantially different from the model value, 
that could produce a systematic error in the inference of 
the global parameters. 

Use of this new model of the Galactic magnetic field, 
which reproduces most large-scale features seen in the ro- 
tation measure and polarized synchrotron skies, should 
allow significant improvements in a number of related 
analyses. With a trustworthy model of the GMF, rota- 
tion measure data can be added to previous constraints 
on the Galactic distribution of thermal electrons ( |Cordes| 
fc Lazio| [2QQ2| . The effects of spatially varying and 
anisotropic diffusion due to the large-scale regular and 
striated GMF can now be included in the determination 
of the Galactic distri bution of cosmic r ays using a code 
such as GALPROP (|Strong et al.||2QQ7[ ). On account of 
the out-of-plane GMF, this may have a significant impact 
on the predicted spectrum and distribution of Galactic 
cosmic rays. Indeed, we have preliminary evidence that 
the typical density of cosmic ray electrons is greater and 
has a larger scale height than predicted by the currently 
standard GALPROP analysis. Finally, a reliable model 
of the large-scale Galactic magnetic field will allow the 
arrival directions of ultrahigh energy cosmic rays to be 
corrected for deflection in the large scale magnetic field, 
for a given charge assignment. 
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APPENDIX 

NOTES ON MODELING ATTEMPTS 

Before arriving at the Galactic magnetic field model described in this paper considerable efforts were made to 
develop and test alternative models. In this appendix we briefly describe the most important of these attempts. We 
also describe a way to implement a more general striated field model, where separate model components (e.g., disk 
and halo) can have different degrees of striation. 



Rejected model features 



i) | Van Eck et al. (2011) presented an extension of the Brown et al. (2007) disk model, based on additional RM data 
in the disk. We im plemented this model and found that it did not improve the y 2 of t he final fit compared to the 
|Brown et al. (2007) model. We chose to base our disk model on the Brown et al. ( 2007[ ) model because it is simpler 
and could be easily modified to conserve magnetic flux. 

ii) Our toroidal halo components initially had the freedom to reverse direction in the inner part of the Galaxy (cf. 
JFWE09). With the inclusion of striated fields, and the subtraction of the local HI bubble described in £ 3.1.1| this 
model feature is no longer necessary in order to explain the observed data; it is sufficient that the outer region has a 
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negligible halo field. 

iii) The field strength in the toroidal halo was given a radial exponential fall-off, but the optimized scale length was 
much larger than the size of the Galaxy. That is, given the radial extent of the electron densities, the observables 
are not sensitive to the outer limits of the GMF. We thus removed this parameter and implemented a constant field 
strength up to radius r n (r s in the south). 

iv) Several different axisymmetric and divergenceless out-of-plane field configurations were tested and optimized, 
including dipole-like models. The final field model was chosen because it gives the best fit to the data of the model- forms 
we considered. 

v) We considered a halo field consisting of axisymmetric or bisymmetric spirals. We let all relevant parameters in 
the spiral fields be free, including the relative orientation in the north/south (i.e., the north and south fields were 
allowed to be completely aligned to completely disaligned). The simpler, toroidal halo field described in the text gave 
a better fit to data, however. 



vi) Striated fields where the level of striation differs for the disk, toroidal halo, and X-field (see SA.2) were also 
considered. No appreciable improvement of the model fit was found. Future work will consider further generalizations, 
such as a purely vertical striated field components (physically motivated by Galactic winds, lifting and stretching field 
lines from the disk) etc. 

Generalized implementation of striated fields 

The simplistic implementation of striated fields done in this paper - by adding a multiplicative factor in the calcula- 
tion of the Stokes' parameters - can be generalized such that different magnetic components (e.g., disk, toroidal halo, 
and X-field) can have different amounts of striation. This i mpl ementation breaks the degeneracy between striated 
fields and a rescaling of the relativistic electron density (see j |5.2|) . 

The most straightforward implementation is simply to include the actual random striated fields explicitly when 
calculating the observables. However, this will often be computationally prohibitive as many realizations are necessary 
to get a reliable mean value. In addition, the stochastic nature of the calculated observables can make the interpretation 
of the MCMC difficult. Instead, we developed a non-stochastic approach that works when the number of different 
striated field components is low. 

As an example we consider the case where the disk, toroidal halo, and X-field components all have a separate 
striated field aligned locally with its regular field. Including any kind of random field makes the calculated observable 
stochastic. As long as only the ensemble average of observables is needed, we can ignore striated fields completely in 
the calculation of rotation measures, since the contribution of random fields to RMs is on average zero. 

To calculate the contribution to the Stokes parameters for a given volumetric cell, let the local magnetic field be 
B regji ± KiB reg ^, where i labels the magnetic field components (disk, toroidal halo, X-field) and the second term 
is the local striated field, which is either parallel or anti-parallel with the regular field. The relative strength of the 
striated field is set by the value K{. In our example, i can take three different values, corresponding to the disk, halo, 
and X-field, and there are thus 2 3 = 8 possible configurations of the local magnetic field. Since we are only interested 
in the average of many realizations, and Stokes parameters are additive, we calculate I^Q^U for a given line-of-sight 
8 times, once for each possible choice implied by adding/subtracting the striated term. We can then simply take the 
mean of the 8 different predicted Stokes parameters. This should correspond to the mean of a very large number of 
stochastic realizations of the field. 
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